import numpy as np
import matplotlib.pyplot as plt


def mackeyglass_eq(x_t, x_t_minus_tau, a, b):
    # Mackey-Glass 方程的右侧表达式
    return a * x_t_minus_tau / (1 + x_t_minus_tau ** 10) - b * x_t


def mackeyglass_rk4(x_t, x_t_minus_tau, deltat, a, b):
    # 使用四阶 Runge-Kutta 方法计算下一个时间步的 x 值
    k1 = deltat * mackeyglass_eq(x_t, x_t_minus_tau, a, b)
    k2 = deltat * mackeyglass_eq(x_t + 0.5 * k1, x_t_minus_tau, a, b)
    k3 = deltat * mackeyglass_eq(x_t + 0.5 * k2, x_t_minus_tau, a, b)
    k4 = deltat * mackeyglass_eq(x_t + k3, x_t_minus_tau, a, b)
    x_t_plus_deltat = x_t + (k1 + 2 * k2 + 2 * k3 + k4) / 6
    return x_t_plus_deltat

def mackeyglass_rk2(x_t, x_t_minus_tau, deltat, a, b):
    # 使用二阶 Runge-Kutta 方法计算下一个时间步的 x 值
    k1 = deltat * mackeyglass_eq(x_t, x_t_minus_tau, a, b)
    k2 = deltat * mackeyglass_eq(x_t + k1, x_t_minus_tau, a, b)
    x_t_plus_deltat = x_t + (k1+ k2)/2
    return x_t_plus_deltat

# 参数设置
a = 0.2  # 方程中的参数 a
b = 0.1  # 方程中的参数 b
tau = 17  # 时滞常数 tau
x0 = 1.2  # 初始条件 x(0) = x0
deltat = 0.1  # 时间步长（与积分步长一致）
sample_n = 60000  # 样本总数，不包括初始条件
interval = 1  # 每隔 'interval' 个时间步输出一次结果
dimension = 1

time = 0.0
index = 0
history_length = int(tau / deltat)
x_history = np.zeros(history_length)  # 假设在 -tau <= t < 0 时，x(t) = 0
x_t = x0

X = np.zeros(int(sample_n/(1/deltat)) + 1)  # 存储所有生成的 x 样本
T = np.zeros(int(sample_n/(1/deltat)) + 1)  # 存储对应的时间样本
j = 1

for i in range(sample_n + 10 ):

    if i%(1/deltat) ==0:
        T[int(i/10)] = time
        X[int(i/10)] = x_t
    if (i % interval == 0):
        print('{:4d} {:f}'.format(i // interval, x_t))
    if tau == 0:
        x_t_minus_tau = 0.0
    else:
        x_t_minus_tau = x_history[index]

    # 计算下一个时间步的 x 值
    x_t_plus_deltat = mackeyglass_rk2(x_t, x_t_minus_tau, deltat, a, b)

    if tau != 0:
        x_history[index] = x_t_plus_deltat
        index = (index + 1) % history_length  # 更新索引，确保循环
    time += deltat

    x_t = x_t_plus_deltat

Y = X.reshape(dimension, -1, order='F')
#保存数据
np.save('Mul_MG_Data1_6000.npy', Y)
print("文件已保存！")

# 绘制结果
plt.figure(figsize=(16, 4))
plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文乱码
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
plt.plot(T[800:1173], X[800:1173])
#plt.xlim([0, T[-1]])
plt.xlabel('t')
plt.ylabel('x(t)')
plt.title('Mackey-Glass 时间序列 (tau={})'.format(tau))
plt.show()
#X = X.reshape((1, -1)).T
